#include <admodel.h>

  extern "C"  {
    void ad_boundf(int i);
  }
#include <model.htp>

model_data::model_data(int argc,char * argv[]) : ad_comm(argc,argv)
{
  Nobs_tot.allocate("Nobs_tot");
  Kw.allocate(1,8,"Kw");
  Data.allocate(1,Nobs_tot,1,6,"Data");
  Kd.allocate(1,Nobs_tot);
  Chl.allocate(1,Nobs_tot);
  Acdm.allocate(1,Nobs_tot);
  Bbp.allocate(1,Nobs_tot);
  lambda.allocate(1,Nobs_tot);
  convert_lambda.allocate(1,Nobs_tot);
}

model_parameters::model_parameters(int sz,int argc,char * argv[]) : 
 model_data(argc,argv) , function_minimizer(sz)
{
  initializationfunction();
  logSigma.allocate(1,8,"logSigma");
  sigma.allocate(1,8,"sigma");
  #ifndef NO_AD_INITIALIZE
    sigma.initialize();
  #endif
  s0.allocate("s0");
  s1.allocate(-1,"s1");
  s2.allocate(-1,"s2");
  loga1.allocate(1,8,"loga1");
  a2.allocate(1,8,0,1.01,"a2");
  eta.allocate("eta");
  pred_Kd.allocate(1,Nobs_tot,"pred_Kd");
  #ifndef NO_AD_INITIALIZE
    pred_Kd.initialize();
  #endif
  a11.allocate(1,8,"a11");
  f.allocate("f");
}

void model_parameters::preliminary_calculations(void)
{

  admaster_slave_variable_interface(*this);
  Kd=column(Data,1);
  Chl=column(Data,2);
  Acdm=column(Data,3);
  Bbp=column(Data,4);
  lambda=column(Data,5);  
  convert_lambda=column(Data,6);
}

void model_parameters::userfunction(void)
{
  f=0.0;
  sigma=exp(logSigma);
  a11=exp(loga1);
  for (int i=1;i<=Nobs_tot;i++)
  {
  pred_Kd(i)=Kw(convert_lambda(i))+Acdm(i)*mfexp(-(s0+s1*Acdm(i)+s2*pow(Acdm(i),2))*(lambda(i)-443.0))+a11(convert_lambda(i))*pow(Chl(i),a2(convert_lambda(i)))+Bbp(i)*pow((lambda(i)/443.0),eta);
  f+=0.5*(Nobs_tot*log(2*M_PI*square(sigma(convert_lambda(i))))+square((log(pred_Kd(i))-log(Kd(i)))/sigma(convert_lambda(i))));
  }
}

void model_parameters::set_runtime(void)
{
  dvector temp1("{200000}");
  maximum_function_evaluations.allocate(temp1.indexmin(),temp1.indexmax());
  maximum_function_evaluations=temp1;
  dvector temp("{1.e-12}");
  convergence_criteria.allocate(temp.indexmin(),temp.indexmax());
  convergence_criteria=temp;
}

model_data::~model_data()
{}

model_parameters::~model_parameters()
{}

void model_parameters::report(void){}

void model_parameters::final_calcs(void){}

#ifdef _BORLANDC_
  extern unsigned _stklen=10000U;
#endif


#ifdef __ZTC__
  extern unsigned int _stack=10000U;
#endif

  long int arrmblsize=0;

int main(int argc,char * argv[])
{
    ad_set_new_handler();
  ad_exit=&ad_boundf;
  arrmblsize = 100000000;
  gradient_structure::set_MAX_NVAR_OFFSET(30000);
    gradient_structure::set_NO_DERIVATIVES();
    gradient_structure::set_YES_SAVE_VARIABLES_VALUES();
  #if defined(__GNUDOS__) || defined(DOS386) || defined(__DPMI32__)  || \
     defined(__MSVC32__)
      if (!arrmblsize) arrmblsize=150000;
  #else
      if (!arrmblsize) arrmblsize=25000;
  #endif
    model_parameters mp(arrmblsize,argc,argv);
    mp.iprint=10;
    mp.preliminary_calculations();
    mp.computations(argc,argv);
    return 0;
}

extern "C"  {
  void ad_boundf(int i)
  {
    /* so we can stop here */
    exit(i);
  }
}
